Non-linear Evolution of f(R) Cosmologies III: Halo Statistics 
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The statistical properties of dark matter halos, the building blocks of cosmological observables 
associated with structure in the universe, offer many opportunities to test models for cosmic accel- 
eration, especially those that seek to modify gravitational forces. We study the abundance, bias 
and profiles of halos in cosmological simulations for one such model: the modified action f(R) the- 
ory. The effects of f(R) modified gravity can be separated into a large- and small-field limit. In 
the large field limit, which is accessible to current observations, enhanced gravitational forces raise 
the abundance of rare massive halos and decrease their bias but leave their (lensing) mass profiles 
largely unchanged. This regime is well described by scaling relations based on a modification of 
spherical collapse calculations. In the small field limit, the enhancement of the gravitational force 
is suppressed inside halos and the effects on halo properties are substantially reduced for the most 
massive halos. Nonetheless, the scaling relations still retain limited applicability for the purpose of 
establishing conservative upper limits on the modification to gravity. 
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I. INTRODUCTION 

In the so-called f(R) class of models (see [1, 2] and 
references therein) cosmic acceleration arises not from an 
exotic form of energy with negative pressure but from a 
modification of gravity that replaces the Einstein- Hilbcrt 
action by a function of the Ricci or curvature scalar R 
[3-5]. 

Cosmological simulations are crucial for exposing the 
phenomenology of f(R) models. In order to satisfy local 
tests of gravity, f(R) models exhibit a non-linear pro- 
cess, called the chameleon mechanism, to suppress force 
modifications in the deep potential wells of cosmologi- 
cal structure [6-10]. Upcoming tests of cosmic acceler- 
ation from gravitational lensing, galaxy and cluster sur- 
veys have most of their statistical weight in the weakly to 
fully non-linear regime. Stringent constraints on modi- 
fied gravity can be expected from current and future sur- 
veys, once the impact on observables in the non-linear 
regime is understood. 

In the previous papers in this series, we have estab- 
lished the methodology for cosmological f(R) simulations 
[11] and conducted a suite of simulations that uncover 
the chameleon mechanism and its effect on the matter 
power spectrum [12]. In this paper, we continue our ex- 
ploration of the non-linear aspects of the f(R) model by 
examining the properties of the basic building blocks of 
cosmological structure: dark matter halos. Specifically, 
we quantify their abundance, i.e. the halo mass func- 
tion, clustering properties, i.e. the linear bias, and their 
density profiles, to see how each are modified from the 



standard cosmological constant, cold dark matter model 
ACDM. 

We begin in §11 with a brief review of the important 
properties of f(R) models and a discussion of the simula- 
tion and analysis methodology. We present our results on 
halo statistics in §111 and discuss them in §IV. Through- 
out we place a special emphasis on exploring the impact 
of the chameleon mechanism and highlighting differences 
between the simulations and conventional scaling rela- 
tions based on linear theory and ACDM. These differ- 
ences expose crucial distinctions that must be considered 
when observationally testing modified gravity theories. 



II. METHODS 

We begin in §11 A by briefly reviewing the basic prop- 
erties of the f(R) model that are important for under- 
standing the cosmological simulations described in §11 B. 
We refer the reader to [12] for a more detailed treatment. 
Finally in §11 C, we discuss the methods used in identi- 
fying the halos and measuring their abundance, bias and 
profiles. 



A. f(R) Gravity 

The f(R) model generalizes the Einstein-Hilbcrt action 
to include an arbitrary function of the scalar curvature 

R, 
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Here L m is the Lagrangian of the ordinary matter and 
throughout c = H = 1 . Force modifications are associated 
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with an additional scalar degree of freedom fn = df /dR. 
For definiteness, we choose the functional form for f(R) 
given in [10] (with n = 1), but neglect higher correc- 
tions of order |/ro| < which results in the following 
effective f(R): 

f(R) = -16TrGp A -f R0 ^. (2) 

Here we define Rq = R(z = 0) and /ho = fn(Ro), where 
overbars denote the quantities of the background space- 
time. For | /ho | <C 1 the background expansion history 
mimics ACDM with = pa/PctH- 

Variation of Eq. (1) with respect to the metric yields 
the modified Einstein equations. We work in the qua- 
sistatic limit, where time derivatives may be neglected 
compared with spatial derivatives. The trace of the mod- 
ified Einstein equations yields the /h field equation 

V 2 <5/h - y [5R(f R ) - SttG^hJ , (3) 

where coordinates are comoving, Sfn — fn(R) — /h(-R), 
SR. = R—R, Sp m = p m —p m . The time-time component of 
the Einstein equations yields the modified Poisson equa- 
tion 

V 2 * = i^a 2 «S Pm -±8R{f R ). (4) 

Here \1/ is the Newtonian potential or time-time metric 
perturbation 2\f r = 5c/oo/.9oo hi the longitudinal gauge. 
These two equations define a closed system for the New- 
tonian potential given the density field. The matter falls 
in the Newtonian potential as usual and so the modifica- 
tions to gravity are completely contained in the equation 
for 

The field equation (3) is a non-linear Poisson-type 
equation, where the non-linearity is determined by 
5R(fn). If the background field /ho is sufficiently large, 
then field fluctuations are relatively small and this term 
may be linearized as SR m (dR/dfn)]^ SJr. It is 
straightforward to show that the Fourier space solution 
to Eqs. (3) and (4) in this approximation is 

rt « = - fe Kr^) ffi2Wk) - (5) 

with p = {idfR/dR)^ 1 / 2 . Hence, gravitational forces are 
enhanced by a factor of 4/3 on scales below p^ 1 , the 
Compton wavelength of the field. We call this regime 
the large field limit. 

Eqs. (3) and (4) in the large field limit imply that the 
field fluctuations are of order the gravitational potential 
\$fit\ ~ 1^1- Therefore if the background field is of order 
the typical gravitational potentials of cosmological struc- 
ture j^l < 10 -5 or smaller, field fluctuations become of 
order unity and SR » (dR/dfji)Sfji which causes the 
Compton wavelength to shrink [10]. We call this the 
small field limit. The large and small field limits are sep- 
arated by a value of the background field of |/ro| ~ 10~ 5 . 



In the small field limit, the field equation (3) then 
requires SR m 8TrG5p m which drives the Poisson equa- 
tion (4) back to its usual form. This is the so-called 
chameleon mechanism which occurs when the back- 
ground field is small compared with the depth of the 
gravitational potential. Hence force law deviations are 
suppressed in the deepest gravitational potentials, i.e. in- 
side the high overdensities of collapsed dark matter halos. 

It is important to note that due to the modified Poisson 
equation (4) for the dynamical potential, the masses dealt 
with in this paper correspond obscrvationally to gravita- 
tional lensing masses, and not to dynamical masses (see 
Appendix A). 



B. Simulations 

To solve the system of equations defined by the modi- 
fied Poisson equation (4) and the /h field equation (3) in 
the context of cosmological structure formation, we em- 
ploy the methodology described in [11] and implemented 
in [12]. Briefly, the field equation for /h is solved on a 
regular grid using relaxation techniques and multigrid it- 
eration [13, 14]. The potential \& is computed from the 
density and /h fields using the fast Fourier transform 
method. The dark matter particles are then moved ac- 
cording to the gradient of the computed potential, — V*, 
using a second order accurate leap-frog integrator. 

We choose a range of background field values |/ho| = 
10~ 6 — 10~ 4 to expose the impact of the chameleon 
mechanism. Since cosmological potentials range from 
10 -6 — 10~ 5 , we expect the chameleon mechanism to be 
operative in the small field limit of this range but ab- 
sent in the large field limit. We also include |/ro| = 
which is equivalent to ACDM. Note that the background 
expansion history for all runs are indistinguishable from 
ACDM to 0(/ho)- More specifically, we take a flat back- 
ground cosmology defined by Qa = 0.76, SI/, = 0.04181, 
Ho = 73 km/s/Mpc and initial power in curvature fluc- 
tuations A s = (4.73 x 10~ 5 ) 2 at k = 0.05Mpc _1 with a 
tilt of n s = 0.958. 

To more directly assess the impact of the chameleon 
mechanism, we also carry out linearized /r simulations in 
which the gravitational potential, ^ , is evaluated accord- 
ing to Eq. (5). In the linearized treatment, the Compton 
wavelength is assumed to be fixed by the background field 
and thus chameleon effects are not present. Therefore, 
the difference between the full /h simulations and the lin- 
earized /h simulations are wholly due to the chameleon 
effects. We will call these runs the "no-chamclcon" sim- 
ulations. 

Tabic I lists the properties of the simulations used in 
the analysis below. All simulations possess 512 grid cells 
in each direction and N p = 256 3 particles. 
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TABLE I: Simulation type and number of runs per box size. 
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C. Halo Properties 

We identify halos and measure their masses in sim- 
ulations with a spherical overdensity algorithm similar 
to [15]. We use cloud- in-cell interpolation to assign the 
particles to the grid. Starting at the highest overdensity 
grid point, we then count the particles within a grow- 
ing sphere centered on the center of mass, until the de- 
sired overdensity with respect to the mean matter density 
A t h = Pm/Pm is reached. Here, we take A t h = 300 for 
defimtencss. The mass A/300 of the halo is then defined 
by the mass of all particles enclosed within this radius 
f 300 ■ We move onto the next highest density grid cell and 
repeat the procedure until all halos have been identified. 
We implicitly take M = M300 below unless otherwise 
specified. 

In our final results we only keep halos with at least 
N m - m dark matter particles, and since our simulations 
are not of high-resolution, we conservatively take N m i n = 
800. We verified that a lower minimum particle number 
of iVmin = 400 provides results consistent with statisti- 
cal uncertainties for all our quoted halo properties. The 
corresponding minimum masses of halos are listed in Ta- 
ble I. 

For each simulation run, we determine the halo mass 
function by binning halos in logarithmic mass intervals, 
and dividing by the comoving volume of the simulation 
box. We then combine different runs and box sizes us- 
ing a bootstrap procedure to produce the estimate of the 
mass function and its errors. We weight each box by 
volume and use only those boxes whose minimum halo 
mass is below the mass bin considered. When measur- 
ing differences between ACDM and f(R), we average the 
differences between simulations with the same initial con- 
ditions to reduce the sample variance. 

We compare simulation results to the Shcth-Tormcn 
(ST) prescription [16] given in Appendix B with modifi- 
cations to spherical collapse as detailed in the Appendix 
A. Semi-analytic prescriptions of this type are widely 
used when analyzing data for cosmological constraints 
and so an assessment of their range of validity is of prac- 
tical importance. 

Next we extract the linear halo bias b^(M ) from our 



simulations. For halos of a given logarithmic mass range 
in a box of size Lbox, we ^ rs t obtain the halo bias b(k 7 M) 
by dividing the halo-mass cross spectrum by the matter 
power spectrum for each simulation 



b{k 7 M) 



Phm(k,M) = (<%(k,Af)S(k)) t 
P mm (fc) (5*(k)<J(k)) k 



(6) 



where #h(k, M) is the halo number density contrast 
whereas <5(k) is the matter mass density contrast. The 
average is over the fc-modes in a fc-bin. For each box we 
employ the modes k > fc m ; n = 2fcf un , where fcf un is the 
fundamental mode of the box (see Table I) and thus the 
smallest boxes barely probe the linear regime. For the 
larger mass bins, we probe more of the linear regime but 
are more limited by small statistical samples. Note that 
the definition of bias adopted will differ from alternate 
choices such as (Phh/Pmm) 1 ' 2 or Phh/Phm in the non- 
linear regime where the correlation coefficient between 
halos and matter can differ from unity. 

In order to remove trends from the non-linearity of the 
bias, we fit a linear relation to b(k, M) = ao(M)+ai(M)k 
between fc min and 10fc m i„, where b(k,M) is the com- 
bined measurement from all boxes. The linear halo bias 
in this mass range is then extrapolated as b-^(M) = 
b(k = 0,M) = ao(M). When considering the modi- 
fications in the f(R) simulations, the same bootstrap 
and linear fit procedure is applied but to the quantity 
Ab/b= (6/(h)-6acdm)/&acdm- Again we compare these 
results to the peak-background split predictions based on 
the ST mass function detailed in Appendix B. 

Finally, we stack the halos in each mass interval and 
measure the average density profile and mass correlations 
of the halos. To reduce scatter within the mass bin we 
scale each density profile to its own r3oo before stacking, 
i.e. we measure 



S P (r/r 300 ) 



Ph{r/r 30 o) 



(7) 



The spatial resolution of our particle-mesh simulations is 
limited by the fixed size of grid cells r ce ii (see Tab. I). We 
measure halo profiles down to the grid scale, though we 
expect that profiles have converged only at scales of sev- 
eral grid cells. When the resolution becomes too low, the 
inner profile flattens leading to a misestimation of both 
the mass enclosed at T3oo and the shape of the halo pro- 
files. We therefore use only the highest resolution boxes 
for our comparisons with the f(R) simulations. The max- 
imum radius for each profile is set to 0.4Lbox- 

In order to avoid biases from incompleteness effects, 
we further limit the range of the stacked profile to radii 
where more than 90% of the halos in the mass bin con- 
tribute. We then bootstrap over all halos in the given 
mass range in order to determine the average profile and 
its error. Profiles and the halo-mass correlation function 
results are compared to the Navarro- Frenk- White (NFW) 
profile and halo model respectively (see Appendix B). 



4 



0.01 ET 




10 12 10 13 10 14 10 15 

M 30 o [M /h] 



FIG. 1: The halo mass function as a function of M300 measured 
in ACDM simulations with bootstrap errors on the mean. The 
upper panel combines different box sizes from 64 to 400Mpc/h and 
compares results with the Sheth-Tormen prediction rescalcd from 
M v to Af3oo as described in the text. The lower panel shows the 
relative deviations from this prediction separately for different box 
sizes. 



III. RESULTS 

In this section we present the results obtained from 
N-body simulations of the f(R) models for the halo 
mass function (§111 A), halo bias (§111 B), density profiles 
(§111 C) and matter power spectrum (§111 D). In all cases, 
we compare the simulation results with predictions using 
scaling relations based on spherical collapse calculations, 
the Press-Schechter prescription and findings from simu- 
lations of ACDM. These calculations are detailed in the 
Appendices. 

Since spherical collapse predictions depend on the 
gravitational force modification, we give a range of pre- 
dictions in each case. The extremes are given by collapse 
with standard gravity and with enhanced forces through- 
out. The former follows the ACDM expectation of a lin- 
ear density extrapolated to collapse of 5 C = 1.673 and a 
virial overdensity of A v = 390; the latter modifies these 
parameters to S c = 1.692 and A v = 309 as detailed in 
Appendix A. 

Neither assumption for the nonlinear collapse is com- 
pletely valid given the evolving Compton wavelength and 
the chameleon mechanism. Moreover, the evolution of 
linear density perturbations used as the reference for the 
scaling relations in Eqs. (Bl), (B4), (B8), and (BIO) as- 
sumes in both cases the full linear growth of the f(R) 
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FIG. 2: Relative deviations of the f(R) halo mass functions from 
ACDM, with \f R0 \ = 1CT 4 (top panel), 10~ 5 (middle panel), and 
(lower panel). In each case, blue squares denote the full 
simulations, while red triangles (displaced horizontally for visibil- 
ity) denote the no chameleon simulations. The shaded band shows 
the range of enhancement expected from spherical collapse rescalcd 
from M v to M3oo- 

model through a(M), including the effects of the evolving 
background Compton wavelength but not the chameleon 
mechanism. Thus unmodified spherical collapse parame- 
ters do not equate to unmodified spherical collapse pre- 
dictions. 



A. Mass Function 

In Fig. 1, we show the halo mass function measured 
from our suite of ACDM simulations along with the boot- 
strap errors described in §11 C. For reference, we compare 
the simulations to the Shcth-Tormen (ST) mass function 
of Eq. (Bl). The ST formula gives the mass function in 
terms of the virial mass and we rescale it to M300 as- 
suming an NFW profile (see Appendix B). Our ACDM 
simulations are consistent with the 10-20% level of accu- 
racy expected of the ST formula and internally between 
boxes of differing resolution. 

Next, we compare the f(R) and ACDM simulations. 
Our measurement of the halo mass function itself is lim- 
ited by statistics and to a lesser extent, resolution (see 
Fig. 1). However, we can reduce the impact of both ef- 
fects by considering the relative difference between the 
halo mass functions measured in f(R) and ACDM simu- 
lations with the same initial conditions and resolution. 
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FIG. 3: The halo bias as a function of wavenumber k in ACDM. 
The upper panel combines different box sizes and runs for halos 
with mass M 30 o = 10 13 - lO 13 - 5 ft -1 M0. The black solid line 
indicates a linear fit, whose extrapolation to k = gives &l (dotted 
red line). Error bars denote bootstrap errors on the mean. The 
lower panel shows the relative deviations from the fit separately 
for each box contributing in this mass range. 



FIG. 4: The linear halo bias as a function of M300 extrapolated 
from the ACDM simulations with bootstrap errors on the mean. 
The upper panel combines different box sizes and runs and com- 
pares the result to the Shcth-Tormen prediction rescaling masses 
from M v to M300. The lower panel shows the relative deviations 
from this prediction. 



Fig. 2 shows this relative enhancement of the halo 
mass- function in the f(R) simulations for different values 
of fm, the background field today, combining different 
box sizes as described in section II C. We show results 
for the full simulations as well as the no-chamcleon sim- 
ulations to help highlight the impact of the chameleon 
mechanism. 

For the large field value of |/ro| = 10~ 4 , the num- 
ber of halos increases significantly, especially at the high 
mass end, by up to 50—150% for cluster-sized halos. 
The chameleon effect slightly suppresses the abundance 
in the high mass end. A similar effect occurs for the 
power spectrum [12] and arises due to the appearance of 
the chameleon effect in deep potentials at high redshifts 
where the background field values are smaller. The over- 
all trend is captured by the spherical collapse predictions 
(shaded band in Fig. 2). The upper limit corresponds to 
unmodified forces, whereas the lower limit corresponds 
to enhanced forces during the entirety of the collapse. 
The enhancement of the linear a(M) in f(R) effectively 
makes objects of the same mass less rare and causes the 
increase in the Sheth-Tormen predictions for the expo- 
nentially suppressed high-mass end of the mass function 
(y = S C /(T > 1). Compared to this effect, that of mod- 
ifying spherical collapse parameters is much smaller. It 
mainly arises from the increase in virial mass with re- 



spect to M300 making the same M300 correspond to rarer 
virialized objects. In this large field limit, all but the 
most massive halos are better described by the modified 
collapse parameters. Moreover, for the purposes of estab- 
lishing upper limits on |/ro| using the halo mass function, 
use of this prediction would only err on the conservative 
side. 

When the value of the Jr field becomes comparable 
to the cosmological potential wells, the chameleon effect 
starts to operate. This can be seen in the mass func- 
tion deviations for \fno\ = 10~ 5 and 10~ 6 (see Fig. 2). 
For the smallest field value, the departures from ACDM 
become very small, so that individual high-mass halos 
change only slightly in mass. Due to the limited statis- 
tics in our simulation sample, we are not able to reliably 
estimate the uncertainties on the mass function deviation 
for the highest mass bin in this case. However, the mean 
deviation in this mass bin is consistent with zero. 

The no-chameleon simulations show a behavior of in- 
creasing deviations at high masses similar to the large- 
field case, while the full f(R) simulations deviate signif- 
icantly from this trend, especially at high masses. For 
|/-f?o| = 10~ 6 the excess almost entirely disappears at the 
highest masses leaving a pile-up of halos at intermediate 
masses. As in the power spectrum [12], the chameleon 
mechanism qualitatively changes the predictions for the 
mass function for |/,ro| i5 10 -5 . 
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FIG. 5: Relative deviations in the halo bias, Ab/b = — 
^ACDm)/^ACDM > as a function of wavenumber k between \Jrq\ = 
W~ 4 and ACDM for M 300 = 10 13 - IO^H^Mq. The black 
solid line indicates a linear fit to the bootstrap means and errors 
of the combined boxes, whose extrapolation to k = gives Ab^/b^ 
(dotted red line). 



FIG. 6: Relative deviations in the f(R) linear halo bias from 
ACDM, with |/ro| = 1CT 4 (top panel), 10~ 5 (middle panel) and 
10 — 6 (lower panel). The no chameleon simulations are again dis- 
placed horizontally for better visibility. The shaded bands show 
the range of deviations of halo bias in f(R) expected from spherical 
collapse with the upper limit corresponding to modified spherical 
collapse parameters. 



It is also apparent from Fig. 2 (lower panel) that the 
spherical collapse predictions are less accurate for the 
small field limit. The range of predictions encompasses a 
deficit of high mass halos that is not seen in the simula- 
tions. Since a(M) is calculated from the linear prediction 
at a radius that encloses the mass M at the background 
density, there would be no predicted enhancement of lin- 
ear fluctuations if this radius is larger than the Compton 
scale in the background. This is in spite of the fact that 
in the no-chameleon simulations forces are still enhanced 
once the perturbation collapses to smaller scales. Com- 
bined with the rescaling of the virial mass, this can pro- 
duce a deficit of predicted objects at a fixed overdensity. 
This problem highlights the difficulties in applying scal- 
ing relations between the linear and non-linear regime, 
which were developed for scale-free ACDM type models, 
to modified gravity theories. 

In the case of the full f{R) simulations, the prob- 
lem is partially compensated by the appearance of the 
chameleon mechanism which also reduces the abundance 
of the highest mass objects by eliminating the extra force 
during the collapse. While the full simulation results lie 
within the range of spherical collapse predictions at the 
high mass end, spherical collapse fails to predict the pile 
up of halos at intermediate masses. 

Still, the ST mass function predictions can be used 
to conservatively place upper limits on \fm\ from the 



abundance of halos with M > 10 14 M Q //i. Employing 
the modified collapse prescription for the enhancement or 
zero, whichever is greater, will always underestimate the 
true enhancement in the suite of models we have tested. 
This underestimate becomes a small fraction of the total 
enhancement for \fm\ > 10 -5 . 



B. Halo Bias 

The halo bias computed from Eq. (6) in the ACDM 
simulations is shown in Fig. 3 for halos with masses in the 
range M 300 = 10 13 - 10 135 h^ 1 M Q as an example. The 
points and error bars are bootstrap averages and errors 
of individual bias computations from the various boxes 
and runs. In this case, only boxes with size Lbox = 64 
and 128 /i _1 Mpc have halos in the mass range and con- 
tribute to the bias calculation (see Tab. I). Note that 
due to the limited halo statistics and our small simulation 
sample, the scatter in the errors themselves is significant. 
We have verified that consistent results are also obtained 
with a lowered N m i n = 100 — 400 which increases halo 
statistics allowing the larger, more linear, boxes to be 
used for the bias. In the lower panel of Fig. 3 we show 
the variation of the bias measurements with box size. 
In the regime of mutual applicability, the bias measure- 
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r / F 300 

FIG. 7: Halo density profile, expressed as the fractional over- 
density S p , for M = 10 14 - 10 15 M Q /h measured in the ACDM 
simulations (upper panel). The halo- mass correlation predictions 
(shaded) represent the range with to without (dotted) profile trun- 
cation [Eq. (Bll)] averaged over the same mass bin. The lower 
panel shows the relative deviation and bootstrap errors measured 
in the different boxes from the prediction without truncation. 



ments between boxes are consistent within the statistical 
uncertainties. 

In Fig. 4, we show the linear halo bias in our ACDM 
simulations as a function of halo mass, measured as de- 
scribed in §11 C. We compare these results to the ST 
bias prediction of Eq. (B4). We again remap the virial 
mass M v to M 30 o and plot the prediction for b^(M 300 )- 
The simulation results are consistent within ~ 20% of 
the prediction. 

Whereas the abundance of halos can be significantly 
changed in f(R), their clustering properties arc relatively 
less affected compared with ACDM. In Fig. 5 we show the 
relative difference between the halo bias in f(R) simula- 
tions with \/ro\ = 1CP 4 and ACDM for the same mass 
bin of Fig. 3. For each box and run contribution, we sub- 
tracted the f(R) simulation bias from that of the corre- 
sponding ACDM simulation with same initial conditions 
to form Ab(k,M)/b(k,M). The averages and error dis- 
played are again obtained by bootstrap of the individual 
differences. The same linear fit procedure is applied and 
evaluated at k = to estimate the relative difference in 
the linear bias Ab L {M)/b L {M) = Ab/b(k, M)\ k=0 . 

In Fig. 6 we compare the linear bias from f(R) and 
ACDM simulations, computed as above, and the range 
of predictions from spherical collapse. The bias decreases 
with increasing |/ro| since halos of a fixed mass become 
less rare and thus less highly biased. The chameleon 
effect in the full simulations decreases the difference in 




FIG. 8: Halo density profile S p in the full f(R) (\f R0 \ = 10" 4 , 
colored) and ACDM simulations (black), for different halo masses 
(upper panel). Profiles for 10 13 - 10 14 and 10 14 - 10 15 M Q /h 
have been multiplied by 10 and 100, respectively. The profiles 
of the highest mass halos were obtained from 128 Mpc/h boxes, 
while the lower mass profiles are from 64 Mpc/h boxes. The lower 
panel shows the relative deviation of the f(R) profiles from those 
of ACDM, with bootstrap error bars. 



bias versus the no chameleon simulations as expected. 
As with the mass function, the spherical collapse range 
adequately describes the high mass halos even for the 
small-field chameleon cases due to a fortuitous cancella- 
tion of modeling errors. 



C. Halo Profiles 

The final ingredient in a basic understanding of ha- 
los and cosmological statistics that are built out of them 
is their average profiles. We plot the fractional density 
contrast S p (r/r 30 o) defined in Eq. (7) and measured in 
the ACDM simulations for the largest and hence best 
resolved mass bin in Fig. 7 (upper panel), for different 
box sizes of the ACDM simulations. For reference we 
compare these with the corresponding halo model pre- 
diction (shaded) from the halo-mass correlation function 
of Eq. (Bll), consisting of an NFW profile plus a 2-halo 
term describing the surrounding mass, averaged over the 
same mass bin as the simulations. The range of predic- 
tions shown is bounded from above by a continued NFW 
profile, and bounded from below by an NFW profile trun- 
cated at r v = r39o as used in the halo model description of 
power spectra, Eq. (B8). In the lower panel of Fig. 7, we 
show the same profiles relative to the halo model predic- 
tion with continued profiles. Removing the overall trend 
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with the halo model better reveals the internal consis- 
tency of our simulations. The agreement between the 
smallest box and the larger boxes with coarser resolu- 
tion and smaller particle number is < 20% in case of the 
128Mpc//i boxes, and < 40% for the 256Mpc//i boxes. In 
the following, we show results from the 128Mpc//i boxes 
for the largest mass halos, in order to increase halo statis- 
tics, and from the 64Mpc//i boxes for all other masses. 

Fig. 8, top panel, shows the stacked halo profiles for 
three mass bins, for ACDM and full f(R) simulations 
with | /ro | = 10~ 4 . The lower panel of Fig. 8 shows the 
relative deviation between ACDM and f(R) halo pro- 
files. When scaled to the same overdensity radius, halos 
in ACDM and f(R) apparently have very similar profiles, 
especially in the inner part of the halo. Although a pre- 
cise measurement of the NFW scale radius is not possible 
with our limited resolution, it is apparent that there are 
no dramatic effects of modified gravity on the halo con- 
centration C300 = ^3oo/ r s- Moreover the deviations are 
consistent with zero well within r3oo ■ The same holds for 
the no-chameleon f(R) simulations. 

For the intermediate and larger halo masses, there is 
an enhancement of the halo profile at r/r 3Q0 ~ few, i.e. 
in the transition region between one-halo and two-halo 
contributions. The smallness of the enhancement of ^ m 
can be explained by a partial cancellation between the 
increased linear power spectrum and reduced linear bias 
in f(R) (§ B and § IIIB). However, a quantitative un- 
derstanding of the behavior of the halo-mass correlation 
at these radii is not possible with the simple halo model 
adopted here, as it fails in the transition region between 
one and two-halo terms (see Fig. 7). In the small field 
simulations, the deviations in the halo profiles are too 
small to be measured with our current suite of simula- 
tions. 

Given the relative smallness of the modified gravity 
effects on halo profiles, the main effect of enhanced forces 
in the large field simulations is to change the mass and 
hence the abundance and bias of halos. 



D. Halo Model Power Spectrum 

We can now put the halo properties together and dis- 
cuss statistics that can be interpreted under the halo 
model paradigm outlined in §B. The matter power spec- 
trum P mm is especially interesting in that the enhance- 
ment in the large field f(R) simulations found in [12] was 
not well described by standard linear to nonlinear scal- 
ing relations [22]. Without an adequate description of 
the large field limit, robust upper limits on \fn\, which 
should be available from current observations, are diffi- 
cult to obtain. 

The halo model provides a somewhat more physically 
motivated scaling relation between the linear and non- 
linear power spectra [23]. Specifically we use the same 
range of ST predictions for the mass function and linear 
bias discussed in the previous sections in Eq. (B8). In ad- 




k [h/Mpc] 

FIG. 9: Power spectrum enhancement relative to ACDM for full 
and no-chameleon simulation and different /jjq field strengths. The 
shaded band shows the predictions from the halo model using pa- 
rameters derived from spherical collapse (see text). 



dition, we vary the concentration parameter of the halos, 
using either an unmodified c v (M v ) relation [Eq. (B6)], or 
an unmodified C300 (A/300) = r 30o/ r s- The latter relation 
is motivated by our finding that the inner parts of halo 
profiles are unmodified in f(R) when referred to the same 
overdensity radius (§ IIIC). Converting C300 to the virial 
concentration, we obtain a ~10% higher c v , which in- 
creases the power spectrum enhancement at k > 1/i/Mpc 
through the 1-halo term [Eq. (B8)]. 

The range of halo model predictions is shown in Fig. 9 
for different values of /^o, together with the simulation 
results from [12]. The upper boundary of each shaded 
band corresponds to unmodified spherical collapse pa- 
rameters and unchanged C300, while the lower boundary 
is using the modified spherical collapse parameters, as- 
suming enhanced forces throughout in the f(R) predic- 
tion, and unchanged c v . 

The halo model provides a reasonable approximation 
to the relative deviations in the large field limit out to 
the k ~ 1 — 3 h/Mpc scales that can be resolved by the 
simulations. The modified collapse provide a somewhat 
better and more conservative approximation for the pur- 
poses of establishing upper limits for \fm\ > 10~ 4 . 

The halo model still fails to capture the chameleon sup- 
pression in the small field limit. Its failure is apparent 
even at |/ fl0 | = 10~ 5 for 0.1 < k(h/Mpc) < 1 and is rel- 
atively larger than the error in the mass function, linear 
bias and halo profiles themselves. This range also corre- 
sponds to the regime where the one halo and two halo 
terms are comparable, i.e. where our simple prescription 
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of linear clustering of halos with density profiles trun- 
cated at the virial radius cannot be expected to apply. 

A prescription that seeks to interpolate between mod- 
ified and unmodified force law predictions [23] and a bet- 
ter treatment of the transition regime that includes non- 
linear halo clustering and halo exclusion could potentially 
provide a better description but is beyond the scope of 
this study. 



IV. DISCUSSION 

Dark matter halos are the building blocks of cosmolog- 
ical obscrvablcs associated with structure in the universe. 
Their statistical properties provide many interesting tests 
of cosmic acceleration, especially of those that seek to 
modify gravitational forces. 

Here we have examined the abundance, clustering and 
profiles of dark matter halos in f(R) modified gravity 
models. In these models, gravitational forces are en- 
hanced below the local Compton scale of an extra scalar 
degree of freedom Generically, this extra force leads 
to an enhanced abundance of massive halos and a de- 
crease in the bias of such halos, but relatively little 
change to the density profile or mass correlation around 
halos of fixed mass. The extent of these effects on halo 
statistics depends strongly on whether the background 
scalar field is in the large field (|/ro| ^ 10 -5 ) or small 
field limit (\f R0 \ < 1CT 5 ). 

In the large field limit, forces are modified everywhere 
below the background Compton scale (Ac > 10 Mpc/ft, 
today [12]). The modifications in this regime are rela- 
tively well described by scaling relations for halo statis- 
tics. By modifying spherical collapse parameters to in- 
clude the enhanced forces, we have shown that the mass 
function and linear halo bias can be described well by 
the Sheth-Tormen prescription. The halo-mass correla- 
tion and average density profiles arc little changed from 
ACDM due to a cancellation of effects from the enhanced 
forces and decreased bias. 

Together these provide a description of the enhanced 
matter power spectrum that corresponds to a relatively 
small overestimate of \fm\ by ~ 50% or less. This level 
of accuracy more than suffices for an order of magnitude 
constraint on field values. Moreover, the overestimate de- 
pends only weakly on \Jro\ & n d can largely be corrected. 
In this prescription, concentration uncertainties which 
are unresolved in our simulations should be marginalized. 
Concentration uncertainties also arise from baryonic ef- 
fects in ACDM [24] and marginalization over these leaves 
only the more unique intermediate scale deviations to dis- 
tinguish modifications of gravity [25] . 

In the small field limit, potential wells of dark matter 
halos are comparable to or larger than the background 
fa field, so that the local Compton wavelength decreases 
substantially from the background value. Modifications 
to gravitational forces then decrease in the interior of 
halos by the so-called chameleon mechanism. This de- 



crease has the effect of bringing deviations in all of the 
halo statistics down at the high mass end. At intermedi- 
ate masses, the excess in the halo abundance can actually 
increase further due to a pile up of halos which also sup- 
presses the change in the bias. 

Scaling relations are not as easily modified to include 
the chameleon effect but do still have limited applicabil- 
ity. Due to a fortuitous cancellation of problems associ- 
ated with a small background Compton wavelength and 
the chameleon mechanism, the modified Sheth-Tormen 
mass function can still be used to provide upper limits 
on the field values that err only on the conservative side. 
Likewise the bias description is reasonably accurate for 
intermediate to high mass halos. We caution that this 
fortuitous cancellation docs not apply to all quantities 
that can be built out of halo statistics. For example the 
halo model for the power spectrum ovcrpredicts the en- 
hancement in the weakly non-linear regime. 

To summarize, in the large field limit which encom- 
passes the range that current cosmological observations 
can test, the scaling relations presented here should al- 
ready enable strong tests of the model. However, more 
work in calibrating the effects of f(R) gravity will be re- 
quired when cosmological observations reach the ~10% 
percent level precision required to test the small field 
limit of f(R) modified gravity. 
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APPENDIX A: SPHERICAL COLLAPSE 

In this Appendix, we examine the modifications to 
spherical collapse induced by the enhanced forces of the 
f(R) model and in particular, derive the collapse thresh- 
old 5 C and the virial overdensity A v used in the main 
text. 

We begin with the nonlinear continuity and Eulcr 
equation for a prcssureless fluid of non-relativistic mat- 
ter. When expressed in terms of the gravitational poten- 
tial ^ , these equations are unaltered by the modification 
to gravity that remains a metric theory (e.g. [26]) 



— + -(v-V)v + ffv = — V*, (Al) 
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where 5 = 5p m /p m and spatial coordinates are comoving. 
These can be combined to a second order equation for 5 



d 2 S 



2H 



85 1 <9 2 (1 + 5)v l v' J _ V • (1 + 
dt a 2 dx l dxi a 2 



(A2) 



but require further information about the velocity and 
potential fields to form a closed system. 

The potential is given by the field equation (3) and 
modified Poisson equation (4) in terms of the density 
fluctuation. For the velocity field, we will take an initial 
top hat density perturbation and make the approxima- 
tion that it remains a top hat throughout the evolution. 
This approximation is valid in the limiting cases that the 
Compton radius is either much larger or much smaller 
than the perturbation. 

Given the top hat assumption for the density, the ve- 
locity field in the interior takes the form v = A(i)r to 
have a spatially constant divergence. Its amplitude is 
related to the top hat density perturbation through the 
continuity equation (Al) 



5+-(l + S)A = Q. 
a 



With the relation 

d 2 v l v : > 
dx l dxi 



12A 2 



S 2 



3 (i + sy 



(A3) 



(A4) 



the spherical collapse equation in the top hat approxima- 
tion becomes 



d 2 6 35 4 S 2 



(1+^2 



V^f , (A5) 



which along with Eqs. (3) and (4) (§ II A) complete the 
system. 

We can bring this equation to its more usual form for 
the radius of the top hat by using mass conservation 

M = (47r/3)r 3 p m (l + 5) = const. (A6) 

Therefore the evolution of r and 5 may be related as 



r 



H - 



1 (5 + 25H--- ^ 



3(1 + 5) 



3 1 + 5' 



(A7) 



Combining this relation with the top hat density equation 
(A5), we obtain 



AttG 



[p m + (1 + 3w)p cff ] - ^V 2 $ , (A8) 



where we have expressed the background expansion in 
terms of an effective dark energy contribution. Note that 
this set of equations also applies to any smooth dark en- 
ergy contribution as long as we take 5R = 8irG5p m in 
the Poisson equation. 

For the f(R) system, there are two limiting cases worth 
noting and these both fall into the class of top hat pre- 
serving evolution. In the large field case the Compton 



wavelength is so long that fn ignores the collapse. In this 
case 6R <C %TrG5p m in the interior. In the opposite small 
field case, the Compton wavelength in the background is 
always smaller than the scale of the perturbation. In this 
case 5R = 87rG<5/9 m as in ordinary gravity with smooth 
dark energy The two limits for the top hat equation 
(A5) can be parameterized as 



AirG 



[pm + (1 + 3ui)p e ff] 



AttG 



FSp m (A9) 



with F = 1/3 corresponding to the large field limit and 
F = corresponding to the small field limit or smooth 
dark energy. Note that p m in the first term on the right 
hand side stands for the total matter overdensity, so that 
for F = the top-hat overdensity follows the same equa- 
tion of motion as the background expansion in a smooth 
dark energy model. 

We now specialize this equation for a background ex- 
pansion that is close to ACDM, w = — 1 and p B s = p\. 
Rewriting the time derivatives in term of ' = d/dlna, a 
ACDM background and with y = [r — r^aj ' a^jri 



1 O m a 3 — 2fL 



2 ^\ Of 



(A10) 



2 \ l m a 6 + \ Ia a-i 



with 



5 = 



1 



yai/a+ 1 



(1 + Si) - 1 



(All) 



and Si as the initial density perturbation at a^. Turn 
around occurs when r' = or y' = — a /on and collapse 
occurs when r = or y = —aja%. 

Under the assumption that the initial conditions are 
set during matter domination when <5 ^ 1, linear theory 
says that 5 cx a 1+p where 



5 5/ 24 

P= h -\/H F . 

y 4 4V 25 



(A12) 



The initial conditions are then y — and y' = — 5,(1 + 
p)/3. More generally, the linearization of the continuity 
and Euler equations imply 



5" + 3^5' = ^^F5. 



H 



H 2 



(A13) 



The linear overdensity extrapolated to the collapse epoch 
is then a function of F. For collapse during matter dom- 
ination 5 C = 1.686 for F — as usual and 5 C = 1.706 for 
F = 1/3. In Fig. 10 (lower panel), we show the threshold 
for collapse at z = as a function of f2 m . In particular 
for Q m = 0.24, S c = 1.673 for F = and 5 C = 1.692 for 
F = 1/3. 

To relate spherical collapse with virialized halos, one 
also has to modify the virial theorem for f(R). All the 
steps in the usual derivation of the tensor virial theorem 
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FIG. 10: Spherical collapse parameters. The linear overdcnsity 
extrapolated to the collapse epoch 8 C and the virial overdensity A v 
are modified from the flat ACDM values (F = 0) by the enhanced 
forces during collapse (F = 1/3). 



The different dependence on r changes the virialization 
radius to the extent that Wa is important. Let us define 
the ratio at turnaround 



V 



2p, 



cir 



2fl A 



(l + F) Pm (l + F)n m a- 3 (l + 5) 



(A16) 



The relationship between the virial radius and the 
turnaround radius s = r v /r max can then be obtained 
from inverting 



'/ 



2s- 1 
2s 3 - s 



(A17) 



Note that as 77 — > 0, s — > 1/2 as expected. The effect of 
F is to make the A term less important. 

In Fig. 10 (lower panel), we show the virial overdensity 
for collapse at z = as a function of fi m . In particular, 
for F = the virial overdensity is A v = 390 for collapse 
today and for F = 1/3 it is lowered to A v = 309. 

These modifications also imply that the virial temper- 
ature of halos of a fixed virial mass is proportional to 
1/3 

(1 + F)A V /J and hence increases for F = 1/3. Likewise, 
hydrostatic equilibrium masses or any masses defined dy- 
namically by the velocity dispersion of the matter would 
be larger than lensing masses by a factor of (1 + F). 



from the Boltzmann equation still apply to f(R) since the 
Boltzmann equation (energy momentum conservation in 
the metric) is unchanged (see e.g. [27]). The only change 
is in relating the potential energy to the matter in the 
top hat 



W = - 3 -(l + F)^ 
5 r 



(A14) 



The implications for spherical collapse then remain 
largely unchanged when expressed in terms of the turn 
around radius. During matter domination the scalar 
virial theorem still reads W = —2T and W(r max ) = 
W(r v ) + T(r v ) = W(r v )/2 and so r v = r max /2. The dif- 
ference is in the density evolution in spherical collapse. 
The traditional way of expressing the virial overdensity 
A v is to take the overdensity at r v during the collapse 
p m (r v ) and divide by the average density at the end of 
collapse jO m (r = 0). For collapse in the matter dominated 
limit F — gives the usual A v = 177.6 and F = 1/3 gives 
A v = 143.1. 

These conditions are modified by the acceleration of 
the expansion at low redshifts. Following [28], the met- 
ric effect of A can be considered as providing a poten- 
tial energy per unit mass of wa = — 47rGp ftr 2 /3. In- 
tegrating this up through the top hat we get W\ = 
— (4:TrGp c g/5)Mr 2 . The virial theorem with the com- 
bined potential energy gives 



T 



1, 



-W + W A 



(A15) 



APPENDIX B: SCALING RELATIONS 

In this appendix, we present the scaling relations that 
were used for comparisons with the simulations in sec- 
tion III. For the mass function we use the Sheth-Tormen 
(ST) prescription [16]. Though other, potentially more 
accurate, descriptions for ACDM exist (e.g. [17]), this 
choice enables us to explore the changes expected in the 
f(R) simulations from spherical collapse (see Appendix 
A). We also found a good match to the ST mass function 
in our ACDM simulations (§ III A). 

The ST description for the comoving number density 
of halos per logarithmic interval in the virial mass M v is 
given by 



dn 



dlnMv 



M v v 'dhiM v 



(Bl) 



where the peak threshold v = S c /cr(M v ) and 



vf(v) = A^-av 2 \\ + K 2 )- p ]exp[-ai/ 2 /2] . (B2) 

Here <r(M) is the variance of the linear density field con- 
volved with a top hat of radius r that encloses M = 
47rr 3 p m /3 at the background density 



a 2 (r) 



d 3 k 



\W{kr)\ 2 P^{k), 



(B3) 



where P\^{k) is the linear power spectrum and W is the 
Fourier transform of the top hat window. The normal- 
ization constant A is chosen such that J duf{v) = 1. The 
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parameter values of p = 0.3, a = 0.75, and S c = 1.673 
for the spherical collapse threshold have previously been 
shown to match simulations of ACDM at the 10 — 20% 
level. The virial mass is defined as the mass enclosed 
at the virial radius r v , where A v = 390 in the ACDM 
model. We discuss modifications to these parameters for 
the f(R) model in §111. 

The peak-background split for halos predicts that the 
linear bias of halos should be consistent with the mass 
function. For the ST mass function, the bias is given by 
[16] 



6 L (M V ) = 6(fc = 0,M v ) 
< 2 -l 



1 



av 



2p 



5 C 8 c [l + (a^)P] ■ 
For the halo profiles, we take an NFW form [18] 

PnfwM = 



r/r s {l + r/r s ) 2 ' 



(B4) 



(B5) 



where r s is the scale radius of the halo and the normal- 
ization p s is given by the virial mass M v . Wc parametrize 
r s via the concentration c v = r v /r s given by [19]: 



Cv(M v ,z = 0) = 9 



My 

M* 



-0.13 



(B6) 



where M* is defined via <t(M*) = S c . By assuming an 
NFW form, we can also rescale mass definitions from the 
virial mass M v to M300 as outlined in [20]. We use this 
approach to compare these scaling relation predictions 
to the simulations in §111 since the definition of the virial 
mass varies with cosmological parameters and f(R) mod- 
ifications. For a given halo in ACDM, M300 is slightly 
larger than My. Given that we generally rescale to M 30 o, 
when no specific overdensity is given we implicitly take 
M = M 300 , e.g. 



dn 



n\n M 



din M- 



::oii 



din M v 
din M 



.".mi 



(B7) 



These properties are combined together in the halo 
model which treats cosmological statistics associated 



with structures through the halos that form them (see 
[21] for a review). For example, the matter power spec- 
trum can be decomposed into 1-halo and 2-halo terms, 

Pmm(fc) = I 2 (fc)F L (fc)+P lh (fc), 



P ih (k) 



M 2 

d\nM v n lnMv ^\y{k,M v )\ 2 , (B8) 



where 



IQs) = I dhxM v n mAU ^y{k,M v )b L {M v ) . (B9) 

Pm 



Here, y(k, M) is the Fourier transform of an NFW density 
profile truncated at r v , unless otherwise specified, and 
normalized so that y(k,M) — > 1 as k — > 0. Note that 
with the ST mass function and bias, lim^o I(k) = 1. 

Likewise the halo-mass cross spectrum P^m for an in- 
finitesimally narrow mass bin around M v is given by 



M 

Pbm = b L (M v )I(k)P L (k) + ^y(k,M v ). (B10) 



Note that the Fourier transform of this quantity is the 
halo-mass correlation function, or average mass profile 



61m (r) 



(Ph(r)) 

Pm 



d 3 k 
(2^)3 



Phme 



-ik-x 



d 3 k 

6 L (M V ) / ^-I(k)P L (k)e- 



+ 



(2tt) 3 ~ 
Pnfw(t) 



(Bll) 



For comparison with simulations, wc show the /?nfw term 
with and without the truncation at the virial radius in 
§111 C. Both the overly simplistic treatment of halo pro- 
files and the use of linear halo correlations make our sim- 
ple model inaccurate in the region where the one and two 
halo pieces are comparable. 
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